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Abstract 

We investigate geodesies in non-homogeneous vacuum pp-w&ve solutions and demon- 
strate their chaotic behavior by rigorous analytic and numerical methods. For the partic- 
ular class of solutions considered, distinct "outcomes" (channels to infinity) are identified, 
and it is shown that the boundary between different outcomes has a fractal structure. 
This seems to be the first example of chaos in exact radiative spacetimes. 
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1 Introduction 



Over the past few decades the realization that the behavior of some nonlinear dynamical 
systems is extremely sensitive to initial conditions has changed our view on time evolution in 
physics, biology, and even economics. In the context of general relativity, which is a nonlinear 
dynamical theory par excellence, the first systems for which chaotic behavior of solutions to the 
Einstein equations has been recognized and studied were spatially homogeneous but anisotropic 
cosmological models. The "chaotic cosmology" programme initiated in 1968 by Misner directed 
attention to Bianchi IX (Mixmaster) models. The oscillatory dependence of their cosmological 
scale factors in different spatial directions on time has been studied both analytically and 
numerically in many works (see relevant contributions in [0J, [[J, and references therein). For 
a long time, however, it has been discussed how conclusively and rigorously the existence 
of chaos in these relativistic systems has been demonstrated as gauge invariant measures of 
chaotic behavior are required. Complicated nonlinear effects also occur in systems with coupled 
gravitational and scalar or Yang-Mills fields where the number of degrees of freedom is increased 
compared to purely gravitational systems (see for example @-[ITJ). 



Another important type of problems providing a nonlinear dynamical system in the context 
of general relativity is the study of geodesic motion in a given spacetime. It is well known that 
Newtonian many-body systems are chaotic and it is of primary interest to investigate similar 
situations in Einstein's theory. In particular, the chaotic behavior of geodesies in the relativistic 
analogue of the two fixed-centres problem was examined by Contopoulos fll]], Dettmann et al 
12) and Yurt ever ||13| , who investigated null and timelike geodesies in a spacetime consisting 



of two (fixed) extreme Reissner-Nordstrom black holes. These results were generalized by 



Cornish and Gibbons |14| to the Einstein-Maxwell-dilaton two-centre spacetimes. Bombelli 



and Calzetta |T5[ and Letelier and Vieira Jl6| studied chaotic geodesic motion in perturbed 
Schwarzschild spacetimes. The chaotic behavior of a spinning test particle in Schwarzschild 



spacetime was demonstrated by Suzuki and Maeda [0. Karas and Vokrouhlicky [fLql , Sota 



et al and Vieira and Letelier pCfl , fl21 | studied the behavior of test particles in some 



static axisymmetric spacetimes. Chaotic motion has also been demonstrated in (topologically 
nontrivial) Robertson- Walker spacetimes by Lockhart et al [^2| and Tomaschitz |[23|| . 

In our work |24| we announced that geodesic motion in the well-known class of pp-waves, 
which are vacuum type-iV solutions representing the simplest exact gravitational wave space- 
times, is also chaotic. Here we present a detailed analysis. In the next section we investigate 
the geodesic equations for a non-homogeneous case. The chaotic behavior of timelike, null and 
spacelike geodesies in these spacetimes is then established by analytic and numerical fractal 
methods in sections 3 and 4, respectively. Some concluding remarks are given in section 5. 
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2 Geodesies in pp- waves 



The metric of the widely known class of vacuum pp- waves Pol, plane- fronted gravitational 
waves with parallel rays, can be written in the form 

ds 2 = 2 dCdC - 2 dudv - (f + f) du 2 , (1) 

where f(u, () is an arbitrary function of the retarded time u and the complex coordinate ( 
spanning the plane wave surfaces u = const. The non-trivial curvature tensor components 
are proportional to so that ([I]) represents the Minkowski universe when / is linear in (. 
The case / = d(u)( 2 describes the well-known plane gravitational wave (the "homogeneous" 



pp-wave) which has been thoroughly investigated, see [25] for Refs. This textbook example of 
an exact radiative solution has also been used for the construction of sandwich and impulsive 

waves 
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Here we study motion in non-homogeneous vacuum pp-w&ves. The geodesic equations for 
the metric (ffl) are 

( + ^U 2 = 0, (2) 
ii = U = const , (3) 
v + (/ lC C + f.S) U + l(f + f), u U 2 = 0, (4) 

where dot denotes djdr with r being an affine parameter along the geodesic. In addition, we 
also assume a condition normalizing the tangent to the geodesic such that {7 M {7 M = e where 
e = —1, 0, +1 for timelike, null or spacelike geodesies, respectively. In particular, r is a proper 
time for timelike geodesic observers. This condition can explicitly be written as 



1 

v= u 



(5) 



Here we assume U ^ (for {7 = the geodesic equations can simply be integrated yielding only 
trivial null geodesies ( = (q, u = u , v = vit + vq and spacelike geodesies ( = -7= exp(i0 o )T + Co ; 
it = uq, v = ViT + v , where ( , u , V\, v and O are constants). By differentiating Eq. (|5]) with 
respect to r and using (0) we immediately obtain (|]) which can thus be omitted. It suffices to 
find solutions of (0) since v(t) can then be obtained by integrating Eq. (||) while u = Ut + uq. 

Now we concentrate on the complex equation (0) which has the same form for timelike, 
null and spacelike geodesies. The first integral of (0) for / independent of u is 

a + U 2 1Zef = 2E , (6) 

where E is a real constant. Using @ we can further simplify @ into v{r) = U~ 1 (2E — e/2) t — 
2U J IZe /(C( r )) dr. Moreover, it indicates that the Hamiltonian for the system of equations 
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(|2]) written in two real coordinates x and y introduced by ( = x + iy is 

H=l( P 2 x +p 2 y ) + V(x,y) , (7) 

with a potential V(x, y) = \U 2 lZef. For non-homogeneous pp-w&ve spacetimes given by 
/ ~ £ n > n — 3, 4, • • -, the corresponding polynomial potential is called an "n-saddle" . Its shape 
for n = 3 and n = 5 is shown in Fig. 1. The metric ([!]) of these spacetimes is invariant 
under rotation £ — > ( = ( exp(i27r/n) so that for any geodesic there exist other n — 1 geodesies 
differing only by rotations. There is also a symmetry Q — > C corresponding to y —>■ —y. 

For simplicity, from now on we shall assume the simplest case n = 3. Since any constant 
multiplicative factor of V can be removed by a suitable rescaling of the parameter r, we can 
without loss of generality consider a potential of the form 

V(x, y) = §x 3 - xy 2 , (8) 

called a "monkey saddle". Surprisingly the Hamiltonian (0), (§) is a special case of famous 
Henon-Heiles Hamiltonian |39[ which is known to be a "textbook" example of a chaotic system. 



Here, however, quadratic terms in V are missing. This particular case of the Henon-Heiles 
Hamiltonian has rigorously been investigated by Rod 



3 Analytic demonstration of chaos in pp-w&ves 

Rod described in detail the chaotic [] behavior of orbits in the Hamiltonian system given by (0), 
(H). He concentrated on bounded orbits. These only appear in the positive energy manifolds 
x,y,p x ,p y ) — E > for which the level surfaces of the potential (§) have the form of a 
"monkey saddle" (see Fig. 2). The homogeneity of V guarantees that the orbit structure for 
any two positive values of E is isomorphic modulo a constant scale factor and adjustment of 
time: x — > x = Ax, y — > y = Xy and r — ► f = r/vA results in E — > E = X 3 E. Therefore, 
without loss of generality we can restrict attention to any arbitrary value of E, say E = |. 
Geodesies for all other values of E can simply be obtained by rescaling space coordinates and 
the affine parameter. 

There are three simple geodesies Lj(t), j = 1,2,3, for a given E which will be described 
in detail in the next section. In projection to the (x, ?/)-plane they follow the axes y = and 
y = ±\^3x of the three symmetric channels of the connected region bounded by the level 
curves V(x, y) = E > (for notation see Fig. 2). The boundary consists of three disjoint 
branches given by y = ±^(x 2 — 1/x) with asymptotes x = 0, x = ±\/3y: V\ intersecting the 
x-axis, V 2 above and V 3 below the x-axis, respectively. 

1 He called it "pathological" since his paper jio) was written long before the term "chaos" came into vogue. 
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In order to obtain a description of the structure of bounded orbits Rod first constructed 
three basic periodic orbits Ilj(r) in each of the three channels (Fig. 2) having perpendicular 
intersections with Lj(t). Rod proved that these three periodic orbits are unstable and in fact 
are isolated invariant sets for the flow generated by the equations of motion. He also showed 
the existence of two additional periodic orbits H4 and H5 which follow the same trajectory in 
the (x, j/)-plane, only n 5 (r) = n 4 (— r). 

The analysis of chaotic behavior was then translated in j|0| from the (x, y)-plane projection 
into a topological description of the asymptotic sets for the periodic orbits n 1 ,n 2 ,ri3 in the 
three-dimensional energy manifold H = E. The region in which these bounded orbits occur 
can be decomposed into three cells Rj such that, e.g., .Ri is the compact, connected subset 
of the (x, y)-plane bounded by the lines Di, D 2 and Ei, see Fig. 2 (here the bar denotes 
a projection of the corresponding set from the energy manifold H(x,y,p x ,p y ) = E to the 
(x, |/)-plane). Formally, D\ = {(x, y) G L 2 with y > 0}, D 2 and D 3 are defined analogously 
using the symmetry and Ei denotes the minimum distance line segment between V 2 and V3 
(this intersects the branches V 2 and V3 in points given by y — =br); E 2 , E 3 are again defined 
analogously. Thus, Ej and the two D k are boundaries of the isolating block Rj for the flow 
containing the unstable periodic solution IT, . Each of the three cells Rj contains one such orbit 
and no other bounded orbits; IL- is the only invariant set in Rj, i.e., it is isolated. Note also 
that any orbit that leaves the central region Rq = R\ U R 2 U R3 never enters it again due to 
the structure of the acceleration field in the channels beyond Ej. Such orbits "go to infinity" 
through the j-channel. 

Subsequently, Rod presented a description of orbits asymptotic to the basic periodic orbits 
Uj as t — > ±00. It is appropriate to localize these sets of asymptotic orbits by their intersection 
with the boundaries Dk which are closed topological two-discs (Fig. 3). Any point P in the disc 
represents a point in the phase space section whose spatial coordinates are all points (x, y) G Dk 
while impulses are restricted by the energy condition implying p\ + p 2 y = r 2 = 2[E — V(x, y)\. 
The radius r monotonically increases from to \JlE as the points (x, y) G Dk are taken 
between the intersection of Dk with V and the origin (0,0). Therefore, we can identify any 
point P = {pxiPy) of the disc with an orbit passing through the point (x(r),y(r)) G Dk with 
a given velocity (x,y) = (p x ,p y ). Using this we can indicate in D k important points, areas 
and their boundaries relevant to the chaotic behavior of orbits, see Fig. 3 for D 3 (figures for 
D\ and D 2 can simply be obtained by symmetry). In the figure, any arc a + (j) denotes the 
boundary of points representing orbits that go to Ej (and then to infinity) as r increases; orbits 
given by a + (j) approach IL, as r — ► +00. Similarly a~(j) describes orbits asymptotic to Uj as 
t — > — 00. Points lj denote the intersections that Lj(r) has with the circle of velocities at the 
origin (orbits pointed towards lj continue directly to the j'-channel as r — >• +00). 

The analysis in |4(| showed that the asymptotic sets which are bordered by a ± (i) intersect 
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transversally on Dk thus denning open sets W(m, n) called "windows" and that the intersection 
of the closure of any window with its bounding sets is uncountable. This gives the existence of 
orbits that "connect" the three unstable periodic orbits IT,-: they are homoclinic (asymptotic 
to the same periodic orbit in both time directions) or heteroclinic (asymptotic to two differ- 
ent periodic orbits, one in each time direction). It is the existence of bounded orbits which 
temporarily enter the region Rj, winding past Hj and then emerging for recycling through the 
next region, and the existence of homoclinic and heteroclinic orbits that illustrates complicated 
chaotic structure of the flow in the central region R . 

In order to describe the topology of all possible orbits in phase space it is convenient 
to introduce a set of sequences (finite, infinite, or biinfinite), s = ■ ■ ■ ,Sk, Sfc+i, Sfc+2, • • •, where 
Sk G {1,2, 3}, Sk 7^ Sfc+i. Any sequence then corresponds to a sequence • • • , R Sk , R Sk+1 , R Sk+2 i • ' ' 
of blocks through which the orbit is running in the prescribed order as r goes from — oo to 
+oo. Using this symbolic dynamics, Rod defined thirteen disjoint orbit classes intersecting Rq\ 

1. bounded orbits with s = {sk}Z^ which do not leave Rq 

2. homoclinic orbits with s = si, S2, • • • , Si that come from n si and go through the finite 
sequence to the same II S1 

3. heteroclinic orbits with s — si, s 2 , • • • , s/ that come from n si and go to different U Sf 

4. orbits given by s — s±, S2, • • • that come from II Sl 

5. orbits given by s = ■ ■ ■ , s_2, S-i that go to II S _ 1 

6. orbits that come from II Sl and go through a finite sequence s = s±, S2, ■ • ■ , s/ to S S/ 

7. orbits that come from and go through a finite sequence s — s^, • • • , s_ 2, s_i to n s , 

8. orbits that come from IL, and go to Sj intersecting only the region Rj for j = 1,2,3 

9. orbits that come from Ej and go to Uj intersecting only the region Rj for j = 1, 2, 3 

10. orbits that come from E Sl and go through the sequence s = S\, S2, ■ • ■ 

11. orbits that go through the sequence s — • • ■ , s_ 2) s_i to E s _ 1 

12. unbounded orbits that come from E S1 through a finite sequence s — s±, s 2 , • • • , Sf to E S/ 

13. orbits that come from Ej to the same Ej intersecting only the region Rj for j = 1,2,3 

Clearly, considering the reversibility of time, the classes 5., 7., 9. and 11. are equivalent to 4., 
6., 8. and 10, respectively. Also, the union of orbit classes 1. through 5. together with the 
periodic orbits III, n 2 , il 3 form the invariant set of bounded orbits. Orbits that are unbounded 
in only one time direction are in classes 6. through 11., and those that are unbounded in both 
time directions are in classes 12. and 13. 
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Finally, the chaotic structure of bounded orbits was analysed in |30] by (a topological 
version of) the Smale horseshoe map. With this technique Rod showed that to any biinfmite 
sequence of symbols {1,2,3} there exists an uncountable number of bounded orbits running 
through the blocks Rj in a given sequence. (Similarly, all the orbit classes 4. through 13. each 
contains an uncountable number of orbits.) Also, the flow admits at least a countable number 
of homoclinic and heteroclinic orbits. 

Rod remarked that the results could be refined if the basic orbits Uj were known to be 
hyperbolic so that they would admit stable and unstable asymptotic manifolds. Consequently, 
to each periodic symbol sequence there would correspond a countable collection of periodic 
orbits. The difficulties in proving the hyperbolicity of H, were subsequently overcome (in an 



even more general context) in [^TJ: the "monkey saddle" potential studied above is a special 
case of "Example A" of fy] given by V(x,y) = — PiV) Pi = 0, p 2 = \/3, 

p3 = -V3- 



In fl4"2f , summarizing and generalizing some previous results |[43|| , the Hamiltonian (|7|), (Q) 



(as a particular case of the Henon-Heiles Hamiltonian) was presented as an example of a system 
for which the Smale horseshoe mapping can be explicitly embedded as a subsystem into the 
flow along the nondegenerate homoclinic and heteroclinic orbits to hyperbolic unstable periodic 
orbits. The complex behavior of nearby orbits then implied the nonexistence of global second 
analytic integral. This completed the proof of the chaotic nature of the studied system in the 
sense of a rigorous definition of chaos, cf. |44] . 



4 Numerical demonstration of chaos in ^waves 

The equations of motion resulting from (|7]), (||) have a very simple form 

x = y 2 — x 2 , y = 2xy . (9) 

However, their explicit analytic solutions can only be found in very special cases. Of course, 
there are (unstable) geodesies with E = given by x = = y. Less trivially, letting y = for 
all r, one gets radial geodesies Li{r) through the y = channel (analogous geodesies L 2 (t) 
and L 3 (r) through remaining two channels can simply be obtained by rotation, see Fig. 2). 
Solutions of this type starting at r = from the level curve V\ — E (i.e., x(0) = \^3E, 
x(0) = 0) are given by 



x(t) = ^3E 



(10) 



1 + cn((5r 

where p 2 = 3" 1 / 6 2 ^E and cn is the Jacobian elliptic function with modulus k 2 = (2 + v / 3)/4 
i.e., k = sin(y7j7r). The function x(r) monotonically decreases from x(0) = \^3E across x(rt) = 



12 

to x(t s ) = — oo (where the singularity is located). Interestingly, the proper times and r s 
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are finite. In fact, one can easily calculate that (3t s = 2K(k) m 5.53612629 where K(k) is the 
complete elliptic integral of the first kind, and (3t\, = F(k,(p) ~ 1.84537543 where F(k,(p) is 
the elliptic integral of the first kind with sinyj = 3 -1 / 4 2/(l + y/3). Clearly, r s = 3rj,. 
In particular, for E = the solution with y = is simply 

which describes (for r > r s ) a geodesic emerging from the singularity x = — oo at r = t s and 
approaching the origin x = asymptotically as r — > oo, or (for r < r s ) a geodesic that falls 
from the origin to the singularity. Again, if we set the initial condition x(0) = Xq < then 



t s = J — G/xq so that the time needed to fall from xq to the singularity is finite. 

We also considered to localize the trajectories of the basic periodic orbits IL,. Due to the 
symmetry one can concentrate on Ux only which can be described as a function x(y). Assuming 
E = 3, it must be a solution of a non- linear equation 

l + ^ X f; X \ " + 3xyx' = \{y 2 - x 2 ) , (12) 
1 + x u 

where x' = dx/dy, such that x(—y) = x(y) and x'(0) = 0. Then Hi is given by 

x = a + by 2 + cy A + dy 6 H , (13) 

where b = |a 2 /(a 3 -l), c= ^(2 + 13a 3 )/(a 3 - l) 2 , d = ^a(505a 6 -437a 3 -68)/(a 3 - l) 4 . We 
found the value of a numerically and then calculated the remaining constants using the above 
relations: a = -0.5152, b = -0.1751, c = 0.0107, d = -0.0012. 

The explicit geodesies presented above are very special. In order to get the global picture 
of a motion one has to perform the integration of Eqs. (||) numerically. Typical geodesies in 
the studied spacetime are shown in Figs. 4 and 5. 

In Fig. 4 we present geodesies starting from the branch V2. The curves are orthogonal to 
V2, proceed first downwards and "fan out" from both sides of 111 (and also II3). Such behavior 



illustrates some of the analytic results presented in |4(| and indicates that all Uj are unstable. 

Other geodesies that we obtained by numerical integrations are shown in Fig. 5. Their 
initial conditions are chosen such that the geodesies start at r = from a circle x 2 + y 2 = p 2 in 
the (x, y)-plane. Due to a scaling property of the "monkey saddle" potential (see section 3) we 
can, without loss of generality, assume p — 1 (all other geodesies except those intersecting the 
origin can simply be obtained by rescaling). In Fig. 5a we present geodesies with x = = y 
at r = and in Fig. 5b geodesies starting with non-vanishing (but same) velocities. (Note 
that these conditions are different from the approach adopted in the previous section since the 
geodesies are not on the same "energy manifolds" E = const. However, E is not the energy of 
particles and photons and there is no physical reason to sort all the geodesies according to E 
no matter how useful it proved to be from the mathematical point of view). 
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We observe that each unbounded geodesic escapes to infinity (i.e., falls to the singularity) 
through only one of the three channels in the potential. Choosing non-zero initial velocities 
makes more geodesies prefer one of the channels (cf. Fig. 5a and Fig. 5b) but does not 
significantly change the character of the motion. 

In fact, all unbounded geodesies through the three channels oscillate around the corre- 
sponding "basic" radial geodesies Lj(r) discussed above. Let us assume geodesies through 
the first channel centered by Li(t) given by (|T0| ) (similar results for the second and the third 
channel follow by symmetry). As they approach the singularity at x = — oo, we may assume 
\y\ <C \x\ (this is also justified by our numerical simulations). Then the asymptotic solution to 
Eq. (D) is x(r) « — 6(r s — r)~ 2 and therefore 

y(r) « [A cos ln(r s - r)] + B sin [^f ln(r s - r)]) , (14) 

where A and B are arbitrary constants. As the geodesies approach the singularity, their 
frequency of oscillations around Li(t) grows to infinity while the amplitude of oscillations 
tend to zero. We call this effect a "focusation" . 

The main objective of this section, however, is to establish the chaotic behavior of geodesies 
in pp-waves by numerical means. This may be more illustrative than the formal analysis pre- 
sented in the previous section. Chaos is usually indicated by a sensitive dependence of possible 
outcomes on the choice of initial conditions. The standard approach, called a fractal method, 
was advanced in the papers by Cornish, Dettmann, Frankel, Levin and others (see for example 
[0, [TT|-[|T4J). It starts with a definition of several different outcomes, i.e., types of ends of all 
possible trajectories. Subsequently, a set of initial conditions is evolved numerically until one 
of the outcome states is reached. Chaos is established if the basin boundaries which separate 
initial conditions leading to different outcomes are fractal. As we shall now demonstrate, we 
observe exactly these structures in the studied system. 

It is natural to parametrize the unit circle from which the geodesies in Fig 5a start at r = 

(with vanishing velocities) by x(0) = cos0, y(0) = sin (ft, (p £ [— ir, ir). We have already pointed 

out that all unbounded geodesies approach the singularity at infinite values of x and y through 

only three distinct channels. These represent possible outcomes of our system and we assign 

them symbol j which takes one of the corresponding values, j £ {1,2,3} (thus, for example, 

j = 1 means that the geodesic approaches the infinity at x — — oo, y = through the first 

channel centered by L\{r) as r — > r s > 0). From Fig 5a we observe that the behavior of the 

function j(<f)) depends very sensitively on <p in certain regions. We calculated j(<fr) numerically 

and we display the results in Fig 6 f\. Also, we plot in the same diagram the function r s (<f)) 

which takes the value of the parameter r when the singularity is reached by a given geodesic. 

2 When the values of j were colour coded we obtained nice fractal pictures. Unfortunately, here we could 
present only their black-and-white versions which are not so impressive. Therefore, we chose simply to plot the 
function j ((f>) . 
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The boundaries between the outcomes appear to be fractal. This is confirmed on the 
enlarged detail of the image and the enlarged detail of the detail etc. up to the sixth level. 
In Fig. 7 we show such zooming in of the fractal interval localized around the value 0^0 
(there are two similar fractal intervals around <p ~ I 71 " an d ~ — §7T corresponding to the 
other outcome channels). On each level the structure has the same property, namely that 
between two larger connected sets of geodesies with outcome channels ji and j'2 7^ ji there is 
always a smaller connected set of geodesies with outcome channel j$ such that j$ 7^ j\ and 
J3 7^ J2- Similarly as in [jllj, [141, the structure of initial conditions of these three types of 
orbits resembles three mixed Cantor sets, and this fact is again a manifestation of chaos. 

The above structure of has its counterpart in the fractal structure of r s (</>), see Fig. 6 
and Fig. 7. We observe that the value of this function increases considerably on each discon- 
tinuity of j{4>), i.e., on any fractal basin boundary between the different outcomes (in fact, t s 
is infinite there). Thus, there is an infinite number of peaks, each of which corresponds to an 
orbit asymptotic to one of the periodic orbits Hj as t s — * 00 (these orbits never "decide" on 
a particular outcome, and so never escape to infinity). Also, the value of t s increases in the 
non-chaotic regions of (ft as one zooms into the higher levels of the fractal. This can be under- 
stood physically. Most geodesies fall into the singularity directly. In Fig. 6 they form the three 
largest connected sets of length A</> ~ 1.83 on which j(0) is constant. Three smaller connected 
sets of length A<pi « 0.248 in the first level of the fractal structure correspond to geodesies 
that approach the singularity after one "bounce" on the potential wall Vj in the central region. 
Their values of r s must naturally be larger. (Note that they also contain geodesies Lj(r) given 
by (f) = 0, (f) = |7T and (p = — |7r for j — 1, 2 and 3, respectively; we calculated explicitly their 
times t s = 2K(k) « 5.5361//? ~ 5.1519 which agrees with the values in Fig. 6). In the second 
level of the fractal there are sets of length A0 2 « 8.6 x 10~ 3 with geodesies that reach the 
singularity after two "bounces" having even larger values of t s etc. We found numerically that 
A0 3 « 3.23 x 10~ 4 , A0 4 « 1.18 x 10~ 5 , A0 5 w 4.28 x 10~ 7 , A0 6 w 1.56 x 10~ 8 , etc. 

Using this we can estimate the fractal dimension D of the above structure. We observe 
that A0j/A0j +1 — >• r = 27.4 ± 0.2 for large values of i. Since the "fractal pattern" is doubled 
in each step we get D = In 2/ lnr = 0.209 ± 0.001 so that the dimension is clearly non-integer. 

The motion can also be visualized by the time evolution of a ring of free test particles in 
the (x, y)-plane which are in rest initially. In Fig. 8 we observe that the circle is deformed in a 
complicated way. In fact, the circle forms loops that escape to different outcome channels as r 
grows. They approach the singularity and one can easily check that as the particles forming the 
original circle move in different channels, their relative proper distance grows. Of course, there 
are also isolated particles of the initial ring that will forever remain near the origin approaching 
asymptotically the basic periodic orbits II,. 

Finally, in order to confirm once more the chaotic behavior of geodesies in pp-w&ves we 
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return back to the Rod analysis. We checked his main result concerning the motion in the 
energy manifold E = const, by numerical simulation of the disc D3. This is shown in Fig. 9 
which significantly improves Rod's schematic sketch (cf. Fig. 3). There are two narrow chaotic 
bands with a fractal structure. Each boundary in the band between two different outcomes 
represents some bounded geodesic as r — > 00. Clearly, there is an uncountable number of such 
geodesies. A similar set for r — > — 00 can simply be obtained by a reflection with respect to 
the y-axis. An interesting feature of the disc is that the fractal bands bend as they approach 
the boundary of the disc which represents the geodesies crossing the origin x = and y = 0. 
In Fig. 10 we present some geodesies of this type. Again, the outcome is extremely sensitive 
on initial conditions. It also demonstrates that the higher the level of the fractal, the greater 
is the number of "bounces" that a geodesic undergoes before falling into the singularity. 



5 Final remarks 

In this paper we have demonstrated by invariant analytic and numerical methods the chaotic 
behavior of geodesic motion in non-homogeneous vacuum pp-waye spacetimes. It was estab- 
lished for all types of geodesies: timelike, null and spacelike. This seems to be the first explicit 
demonstration of chaos in exact radiative spacetimes (note that a chaotic interaction of par- 
ticles with particular classes of linearized gravitational waves on given backgrounds has been 
studied in |TB| , P5fl-f4"8fl). In fact, pp-w&ve solutions represent the simplest models of 



exact gravitational waves in general relativity. However, so far most work has concentrated 
on homogeneous pp-w&ves for which the motions are integrable (Eq. ([|) is linear) and hence 
non-chaotic. 

Although we investigated in detail only the spacetime (HD with the function f(() of the 
form / ~ ( 3 corresponding to the "monkey saddle" potential V given by (§), as we indicated 
in |24j the results might be carried out to / ~ £ n with n > 4, i.e., to the general n-saddle 



potential V = ^7Ze( n (cf. Fig. 1 for n = 5). In particular, the decomposition of the central 
region into topological isolating blocks and the existence of basic unstable periodic solutions 
n,- in each channel j = 1, 2, • • • , n is immediate [EO]. These periodic solutions are hyperbolic 



[ f4"T| and the existence of nondegenerate homoclinic and heteroclinic orbits is established 
[ f43| . As an example we present our numerical results for n = 5: Fig. 11a shows unbounded 
geodesies starting from a unit circle from rest and Fig. lib displays corresponding functions 
j((f>) and T s (<f)). There are now five outcome channels and again, the structure is fractal. 

Our work demonstrates that geodesic motion in spacetimes even as simple and well-known 
as pp-waves can be complex. Hopefully, it will initiate investigation of chaos in other exact 
radiative space-times. 
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Figure Captions 



Fig. 1 The shape of the potential V(x, y) — ^ Tie ( n , where ( = x + iy, describing geodesies in 
corresponding non-homogeneous vacuum pp-w&ve spacetimes (here for n = 3 and n — 5). 

Fig. 2 For n — 3, the "monkey saddle" potential V(x, y) = |x 3 — xy 2 plays a crucial role 
in famous Henon-Heiles chaotic system. The level surface V — E consists of disjoint 
branches Vj, j = 1,2,3, defining three channels to infinity, each centered by a special 
geodesic Lj(t). The central region can be decomposed into three cells Rj bounded by 
Dk and Ej such that the only bounded orbit in Rj is an unstable periodic orbit Ilj(r). 
Also, the trajectory of two additional periodic orbits n 5 (r) = n 4 (— r) which we found 
by numerical simulations is drawn by a dashed line. See the text for more details. 

Fig. 3 The boundary D 3 between the isolating blocks R 2 and R 3 is a closed topological two- 
disc. Any point (p x , p y ) in this schematic figure represents an orbit passing with a velocity 
(x,y) = (p x ,Py) through the point (x(r),0) G D 3 such that (for E — |) x — ^Jl — |r 2 
where r = ^jp\ + Py- Arc boundaries of orbits that go to infinity through the j-channel as 
r — > ±oo are denoted by a ± (j). Their intersections form "windows" W(m,n) containing 
homoclinic and heteroclinic orbits asymptotic to and from the basic periodic orbits II. 

Fig. 4 Geodesies starting from the branch V 2 "fan out" from both sides of Tlx and n 3 . 

Fig. 5 Geodesies starting from a unit circle in the (x, y)-plane a) from rest, b) with initial 
velocity x = 0.4, y = 0.8. They escape to infinity (where the singularity is localized) 
through only three channels in the potential and a sensitive dependence of these three 
possible outcomes j G {1, 2, 3} on the choice of initial conditions is observed. 

Fig. 6 Plot of functions j{<p) and T s (<j>) which labels the three possible outcomes and the value 
of r when the singularity is reached by a given geodesic, respectively, on G [— n, ir) 
parametrizing the initial position on a unit circle, x(0) = cos0, y(0) = sin0; x(0) = = 
y(0). Boundaries separating different outcomes are fractal establishing chaos. Each peak 
in t 8 {4>) which coincides with discontinuity in j(0) corresponds to some bounded orbit. 

Fig. 7 The fractal structure of j(4>) and T s ((f>) was confirmed by zooming in the interval around 
the value <fi ~ up to the sixth level. Between any two connected sets representing 
geodesies with outcome channels j± and j 2 ^ ji there is always a smaller connected set 
with j 3 ^ ji and j 3 ^ j 2 . 

Fig. 8 The time evolution of a ring of free test particles in the (x, y)-plane, initially in rest. 
The circle is deformed in a complicated fractal way with different segments moving to 
different outcome channels. Notice, for example, a similarity of the patterns at r = 4.0 
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and r = 7.0. However, at r = 7.0 the lines are in fact doubled and consist of particles 
coming from different parts of the original circle. 

Fig. 9 Exact form of the disc D 3 which we obtained by numerical simulation significantly 
improving Rod's schematic sketch presented in Fig. 3. The chaotic bands in the dies of 
radius « 0.8165 are very narrow and also have fractal structure which is clear from 
the enlarged details. The coding is such that white colour denotes the outcome channel 
j — 1, black corresponds to j — 2 and grey corresponds to j — 3. Each boundary in the 
fractal basin represents a geodesic asymptotic to some IT- as r — > +oo. 

Fig. 10 Typical geodesies starting from the origin x(0) = = y(0) with velocities x(0) = 
yfcos^, y(0) = ^sinip for a) ip e (-0.28, +0.28), b) if> e (0.024,0.028), c) if; e 
(0.0254,0.0256). Again, the outcome is extremely sensitive on initial conditions. 

Fig. 11 Motion in pp-wwes given by the structural function / ~ ( 5 : a) geodesies starting from 
a unit circle from rest approach infinity through five outcome channels, b) corresponding 
functions j(<f>) and T s ((p) also have fractal structure. 
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